### New analysis subsetting votes using Braeuninger et al data
### Produces table A6 in supplementary appendix
### JBS 4 August 2017

rm(list=ls())

library(readstata13)
library(plyr)
library(lme4)
library(stargazer)

# set path here
path <- "~/Dropbox (Personal)/Rebel Summaries/APSR_SKLLO_Repfiles"

# Read in Braeuninger et al data
braundta <- read.dta13(paste(path,"/BraunetalReplication/data/bms_replic_uk.dta", sep= ""))

braundta <- braundta[braundta$term != 2001, c(1,2,3,4,7,8,15,16)]

braundta <- braundta[!duplicated(braundta$ID), ]

braundta$voteid <- paste(braundta$date, braundta$vote_no, sep="_")
braundta$voteid <- paste("X", as.character(braundta$voteid), sep = "")
braundta$voteid <- gsub("-",".",braundta$voteid)


setwd(paste(path,"/RawData",sep=""))

leader <- read.csv("UK HoC Leadership Data 1979-2016.csv")
leader0510 <- leader[,grep("X2005_5",colnames(leader),perl=T):grep("X2010_4",colnames(leader),perl=T)]
leader1015 <- leader[,grep("X2010_5",colnames(leader),perl=T):grep("X2015_4",colnames(leader),perl=T)]

leader0510.mp <- leader0510>0
leader0510 <- leader0510[rowMeans(leader0510.mp)==1  & (leader$Party.Affiliation==1 | leader$Party.Affiliation==2 ) ,]
PubWhipID.05 <- leader$PubWhipID[rowMeans(leader0510.mp)==1 & (leader$Party.Affiliation==1 | leader$Party.Affiliation==2 ) ]

leader1015.mp <- leader1015>0
leader1015 <- leader1015[rowMeans(leader1015.mp)==1 & (leader$Party.Affiliation==1 | leader$Party.Affiliation==2 ) ,]
PubWhipID.10 <- leader$PubWhipID[rowMeans(leader1015.mp)==1 & (leader$Party.Affiliation==1 | leader$Party.Affiliation==2 ) ]


leader0510.bin.ld <- leader0510>=6
leader0510.ld <- rep(0,length(PubWhipID.05))
leader0510.ld[rowMeans(leader0510.bin.ld)==1] <- 1
bin.leader0510 <- as.data.frame(cbind(PubWhipID.05,leader0510.ld))

colnames(bin.leader0510) <- c("PubWhipID","leader.all")
bin.leader0510 <- bin.leader0510[!is.na(bin.leader0510$PubWhipID),]

leader1015.bin.ld <- leader1015>=6
leader1015.ld <- rep(0,length(PubWhipID.10))
leader1015.ld[rowMeans(leader1015.bin.ld)==1] <- 1
bin.leader1015 <- as.data.frame(cbind(PubWhipID.10,leader1015.ld))

colnames(bin.leader1015) <- c("PubWhipID","leader.all")
bin.leader1015 <- bin.leader1015[!is.na(bin.leader1015$PubWhipID),]

dta.0510 <- read.csv(paste(path,"/RawData/2005-2010/rebels.csv", sep=""))
dta.1015 <- read.csv(paste(path,"/RawData/2010-2015/rebels.csv", sep=""))

totdiv.0510 <- ncol(dta.0510)-5
totdiv.1015 <- ncol(dta.1015)-5

## Read in later period free vote data
fv.0510 <- read.csv(paste(path,"/RawData/2005-2010/divisions.csv",sep=""))
fv.1015 <- read.csv(paste(path,"/RawData/2010-2015/divisions.csv",sep=""))

fdiv.0510 <- fv.0510$Division[fv.0510$FreeVote==1|fv.0510$FreeVLab==1|fv.0510$FreeVCon==1]
fdiv.1015 <- fv.1015$Division[fv.1015$FreeVote==1]

fdiv.0510 <- paste("X", as.character(fdiv.0510), sep = "")
fdiv.0510 <- gsub("-",".",fdiv.0510)

fdiv.1015 <- paste("X", as.character(fdiv.1015), sep = "")
fdiv.1015 <- gsub("-",".",fdiv.1015)

# Drop free votes
dta.0510 <- dta.0510[,colnames(dta.0510) %in% fdiv.0510==F]
dta.1015 <- dta.1015[,colnames(dta.1015) %in% fdiv.1015==F]

# Drop opposition sponsored leg and PMBs 

braundta05GOVNOTPMD <- braundta$voteid[braundta$term==2005 & as.numeric(braundta$gov_sponsor)==2 & braundta$privatememb_dy==0]
braundta10GOVNOTPMD <- braundta$voteid[braundta$term==2010 & as.numeric(braundta$gov_sponsor)==2 & braundta$privatememb_dy==0]

dta.0510 <- dta.0510[,colnames(dta.0510) %in% fdiv.0510==F]
dta.1015 <- dta.1015[,colnames(dta.1015) %in% fdiv.1015==F]

### DO SUBSETTING HERE
dta.0510 <- dta.0510[,colnames(dta.0510) %in% braundta05GOVNOTPMD==F]
dta.1015 <- dta.1015[,colnames(dta.1015) %in% braundta10GOVNOTPMD==F]

dta.0510 <- join(dta.0510, bin.leader0510,by="PubWhipID", type="inner") 
dta.1015 <- join(dta.1015, bin.leader1015,by="PubWhipID", type="inner") 

dta.0510 <- dta.0510[,c(1:5,ncol(dta.0510),7:ncol(dta.0510)-1)]
dta.1015 <- dta.1015[,c(1:5,ncol(dta.1015),7:ncol(dta.1015)-1)]



## Calculate rebel data 2005--2010
tot.reb.05 <- rowSums(dta.0510[,7:ncol(dta.0510)],na.rm=T)
tot.absent.05 <- rowSums(is.na(dta.0510[,7:ncol(dta.0510)]))
tot.withparty.05 <- ncol(dta.0510[,7:ncol(dta.0510)])-tot.reb.05-tot.absent.05

tot.reb.10 <- rowSums(dta.1015[,7:ncol(dta.1015)],na.rm=T)
tot.absent.10 <- rowSums(is.na(dta.1015[,7:ncol(dta.1015)]))
tot.withparty.10 <- ncol(dta.1015[,7:ncol(dta.1015)])-tot.reb.10-tot.absent.10

tot.reb.dta.05 <- cbind(dta.0510[,c(1:6)],tot.reb.05,tot.withparty.05)
tot.reb.dta.10 <- cbind(dta.1015[,c(1:6)],tot.reb.10,tot.withparty.10)

merge0510 <- merge(tot.reb.dta.05, tot.reb.dta.10, by="PubWhipID",all=T)
merge0510$mostrebel <- (merge0510$tot.reb.05 + merge0510$tot.reb.10)
merge0510 <- merge0510[,c("PubWhipID","mostrebel")]

tot.reb.dta.05 <- merge(merge0510,tot.reb.dta.05, by = "PubWhipID")
tot.reb.dta.10 <- merge(merge0510,tot.reb.dta.10, by = "PubWhipID")



tot.reb.dta.05$firstterm <- 1
tot.reb.dta.10$firstterm <- 0

tot.reb.dta.05$gov <- 0
tot.reb.dta.05$gov[tot.reb.dta.05$Party=="Lab"] <- 1 

tot.reb.dta.10$gov <- 0
tot.reb.dta.10$gov[tot.reb.dta.10$Party=="Con"] <- 1

median.reb.lab.05 <- median(tot.reb.dta.05$tot.reb.05[tot.reb.dta.05$Party=="Lab"])
median.reb.con.05 <- median(tot.reb.dta.05$tot.reb.05[tot.reb.dta.05$Party=="Con"])

tot.reb.dta.05$reb.med.05 <- 0
tot.reb.dta.05$reb.med.05[tot.reb.dta.05$tot.reb.05 > median.reb.lab.05 & tot.reb.dta.05$Party=="Lab" ] <- 1
tot.reb.dta.05$reb.med.05[tot.reb.dta.05$tot.reb.05 > median.reb.con.05 & tot.reb.dta.05$Party=="Con" ] <- 1

q75.reb.lab.05 <- quantile(tot.reb.dta.05$tot.reb.05[tot.reb.dta.05$Party=="Lab"],0.75)
q75.reb.con.05 <- quantile(tot.reb.dta.05$tot.reb.05[tot.reb.dta.05$Party=="Con"],0.75)

tot.reb.dta.05$reb.q75.05 <- 0
tot.reb.dta.05$reb.q75.05[tot.reb.dta.05$tot.reb.05 > q75.reb.lab.05 & tot.reb.dta.05$Party=="Lab" ] <- 1
tot.reb.dta.05$reb.q75.05[tot.reb.dta.05$tot.reb.05 > q75.reb.con.05 & tot.reb.dta.05$Party=="Con" ] <- 1

median.reb.lab.10 <- median(tot.reb.dta.10$tot.reb.10[tot.reb.dta.10$Party=="Lab"])
median.reb.con.10 <- median(tot.reb.dta.10$tot.reb.10[tot.reb.dta.10$Party=="Con"])

tot.reb.dta.10$reb.med.10 <- 0
tot.reb.dta.10$reb.med.10[tot.reb.dta.10$tot.reb.10 > median.reb.lab.10 & tot.reb.dta.10$Party=="Lab" ] <- 1
tot.reb.dta.10$reb.med.10[tot.reb.dta.10$tot.reb.10 > median.reb.con.10 & tot.reb.dta.10$Party=="Con" ] <- 1

q75.reb.lab.10 <- quantile(tot.reb.dta.10$tot.reb.10[tot.reb.dta.10$Party=="Lab"],0.75)
q75.reb.con.10 <- quantile(tot.reb.dta.10$tot.reb.10[tot.reb.dta.10$Party=="Con"],0.75)

tot.reb.dta.10$reb.q75.10 <- 0
tot.reb.dta.10$reb.q75.10[tot.reb.dta.10$tot.reb.10 > q75.reb.lab.10 & tot.reb.dta.10$Party=="Lab" ] <- 1
tot.reb.dta.10$reb.q75.10[tot.reb.dta.10$tot.reb.10 > q75.reb.con.10 & tot.reb.dta.10$Party=="Con" ] <- 1



names(tot.reb.dta.05) <- c("PubWhipID","mostrebel","Name","MPName","Constituency","Party","leader","tot.reb","tot.withparty","firstterm","gov","med","q75")    

names(tot.reb.dta.10) <- c("PubWhipID","mostrebel","Name","MPName","Constituency","Party","leader","tot.reb","tot.withparty","firstterm","gov","med","q75")    



maj05 <-read.csv(paste(path,"/RawData/2005-2010/votes.csv",sep=""))
maj05 <- maj05[,c("PubWhipID","Majority")] 
maj10 <-read.csv(paste(path,"/RawData/2010-2015/votes.csv",sep=""))
maj10 <- maj10[,c("PubWhipID","Majority")] 

tot.reb.dta.05 <- merge(tot.reb.dta.05,maj05, by="PubWhipID")
tot.reb.dta.10 <- merge(tot.reb.dta.10,maj10, by="PubWhipID")


tot.reb.dta.0510 <- as.data.frame(rbind(tot.reb.dta.05, tot.reb.dta.10))

tot.reb.dta.0510$Majority <- tot.reb.dta.0510$Majority/10000

startyear <- read.csv(paste(path,"/RawData/startdates.csv",sep=""))

startyear$start_year10 <- (2010 - startyear$start_year)/(2010 - min(startyear$start_year))

start10 <- startyear[,c("PubWhipID","start_year10")]

colnames(start10) <- c("PubWhipID","start_year")

tot.reb.dta.0510 <- join(tot.reb.dta.0510,start10, by="PubWhipID", type="left")

#save(braun.dta.0510,file="~/Dropbox (Personal)/Rebel Summaries/Cleaned Data and Scripts for Analysis/calcdata/braunrebeldta0510.rda")

mod1 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+start_year+leader+(1|PubWhipID), family=binomial(link=logit), data= tot.reb.dta.0510[tot.reb.dta.0510$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod2 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+start_year+leader+(1|PubWhipID), family=binomial(link=logit), data= tot.reb.dta.0510[tot.reb.dta.0510$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

summary(mod1)
summary(mod2)


groupsdta <- read.csv(paste(path,"/RawData/groupmembership.csv",sep=""))

dta.w.groups0510 <- join(tot.reb.dta.0510,groupsdta,by="PubWhipID", type="left")
dta.w.groups0510$LCG[is.na(dta.w.groups0510$LCG)] <- 0
dta.w.groups0510$CORNER[is.na(dta.w.groups0510$CORNER)] <- 0
dta.w.groups0510$MONDAY[is.na(dta.w.groups0510$MONDAY)] <- 0
dta.w.groups0510$NTB[is.na(dta.w.groups0510$NTB)] <- 0
dta.w.groups0510$CONGROUPANY[is.na(dta.w.groups0510$CONGROUPANY)] <- 0
dta.w.groups0510$CONGROUPSUM[is.na(dta.w.groups0510$CONGROUPSUM)] <- 0

dta.w.groups0510$group <- 0
dta.w.groups0510$group[dta.w.groups0510$CONGROUPANY == 1 | dta.w.groups0510$LCG == 1] <- 1


mod3 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups0510[dta.w.groups0510$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod4 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups0510[dta.w.groups0510$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

summary(mod3)
summary(mod4)

########################

load(paste(path,"/CalcData/positions_jsnew.Rda",sep=""))

ws.pos <- results1[results1$firstperiod==0,]

# Keep only Labour for now
reblab05 <- dta.w.groups0510[dta.w.groups0510$Party=="Lab",]

# Load rebel data and create name var for merging

fname <- as.character(reblab05$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(reblab05$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

reblab05$fullname <- paste(fname, lname, sep=" ")

## Fix names
## Changed Alan Williams. Only one Alan Williams in 2nd period

ws.pos$names[ws.pos$names=="Alan WilliamsSwansea"] <- "Alan Williams"
ws.pos$names[ws.pos$names=="Huw IrrancaDavies"] <- "Huw Irranca-Davies"
ws.pos$names[ws.pos$names=="Robert MarshallAndrews"] <- "Robert Marshall-Andrews"

mergedta05 <- merge(reblab05 , ws.pos, by.x = "fullname", by.y = "names")

mergedta05$Score <- -1*mergedta05$Score

mod5 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta05)
summary(mod5)

load(paste(path,"/CalcData/positions_tories.Rda",sep=""))

rebcon05 <- dta.w.groups0510[dta.w.groups0510$Party=="Con",]

ws.pos <- results1[results1$firstperiod==0,]

# Load rebel data and create name var for merging

fname <- as.character(rebcon05$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(rebcon05$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

rebcon05$fullname <- paste(fname, lname, sep=" ")


## Fix names
## Changed Alan Williams. Only one Alan Williams in 2nd period

ws.pos$names[ws.pos$names=="David HeathcoatAmory"] <- "David Heathcoat-Amory"
ws.pos$names[ws.pos$names=="Geoffrey CliftonBrown"] <- "Geoffrey Clifton-Brown"
ws.pos$names[ws.pos$names=="Ian LiddellGrainger"] <- "Ian Liddell-Grainger"
ws.pos$names[ws.pos$names=="Ian LiddellGrainger"] <- "Ian Liddell-Grainger"

mergedta05 <- merge(rebcon05 , ws.pos, by.x = "fullname", by.y = "names",all.y=T)


### Models for 05

mod6 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta05)
summary(mod6)

stargazer(mod1,mod2,mod3,mod4,mod5,mod6,type="latex",title="Effect of Government Status on Rebellion: Fractional Logit Models with Random Effects, Government Bills", align=T, column.labels=c("Labour","Conservative","Labour","Conservative"),style="ajps")
